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Calculating the gravitational self force in Schwarzschild spacetime 
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We present a practical method for calculating the local gravitational self-force (often called 
"radiation-reaction force") for a pointlike particle orbiting a Schwarzschild black hole. This is 
an implementation of the method of mode-sum regularization, in which one first calculates the (fi- 
nite) contribution to the force due to each individual multipole mode of the perturbation, and then 
applies a certain regularization procedure to the mode sum. Here we give the values of all the 
"regularization parameters" (RP) required for implementing this regularization procedure, for any 
geodesic orbit in Schwarzschild spacetime. 
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The problem of determining the self force (SF) act- 
ing on a charged particle dates many decades back to 
the classical works by Lorentz jl) and Dirac Q on the 
electron's equation of motion. This force — often referred 
to as "radiation reaction force" — may be attributed to 
the flux of energy-momentum carried away by the elec- 
tromagnetic radiation. Recent years' growing interest in 
sources of gravitational waves has renewed interest in this 
old problem, in a new and exciting context: the gravita- 
tional SF acting on a point mass in curved spacetime. 

This manuscript deals with the motion of a small point- 
like particle of mass fj, around a Schwarzschild black hole 
with mass M ^ jj.. When /i is finite, the particle's 
trajectory deviates from a geodesic of the Schwarzschild 
background geometry, due to the particle's interaction 
with its own gravitational field. This deviation indi- 
cates the presence of a self force. The particle's equa- 
tion of motion thus takes the form iiUa = ^'a*''^ where 
tta is the covariant acceleration (i.e., Oq = Ua-^pu^), and 
F^^^ oc 0(/i^) describes the leading-order gravitational 
SF effect. In general, knowledge of F^^^ is essential 
for modeling the orbital evolution in binary black-hole 
systems with extreme mass ratios. Such astrophysical 
systems — typically, the scenario of a small compact ob- 
ject captured by a supermassive black hole (of the kind 
now believed to reside in the cores of many galaxies, in- 
cluding our own) — are expected to serve as main targets 
for future space-based gravitational wave detectors d). 
The need for precisely predicting the orbital evolution of 
such binary systems — and the waveform of the emitted 
gravitational radiation — strongly motivates the develop- 
ment of practical methods for gravitational SF calcula- 
tions. 

When dealing with point particles, one faces the fun- 
damental issue of regularization: extracting the "correct" 
finite part of the (otherwise divergent) SF. A formalism 



for regularizing the electromagnetic SF in curved space- 
time was developed long ago by DeWitt and Brehme [|| . 
Recently, Mino, Sasaki, and Tanaka (MST) [|| first ob- 
tained a formal expression for the gravitational SF, based 
on a physically well-motivated regularization technique. 
The same formal result was independently derived by 
Quinn and Wald (QW) Q using an axiomatic approach. 
It is useful to write this result as 



lim [F^^^\x)-F^"{x)\ 



(1) 



where Zg is the SF's evaluation point, is a point in the 
neighborhood of Zg , F^"" is the "full" force constructed 
from the metric perturbation (as described below), and 
F^" (the "direct" force) is the "divergent piece" to be 
removed, whose formal construction is described in [^,^ . 
Roughly speaking, the direct piece F^^^ represents the 
instantaneous effect of waves propagating directly along 
the particle's light cone; hence, the SF effect is entirely 
attributed to the so-called "tail" part of the particle's 
field, i.e., waves scattered inside the particle's past light 
cone due to the background's curvature. 

The above regularization prescription by MST is based 
on the physically well-motivated method of matched 
asymptotic expansions. At present, this prescription — 
which also conforms with QW's prescription — is widely 
accepted as the standard regularization scheme for the 
SF 0. As it stands, however, it is difficult to directly 
implement this method in practical calculations, because 
it requires the computation of the Green's function from 
any worldline point to any other future worldline point 
1^ . To allow the practical implementation of Eq. (|l|) in 
strong-field calculations, Barack and Ori (BO) devised a 
method based on a mode decomposition of the pertur- 
bation field (see |9j for the scalar SF and for the 
gravitational SF). In the mode-sum method, the quanti- 
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ties -F^"" and F^" are first formally expanded into mul- 
tipole Z-modes F^^^ and F^". Most benificial is the fact 
that, whereas F^^^ diverges as x — s- zq, its individual 
modes attain finite values even at the particle's location. 
The full mode F^"'' is directly derived (see below) from 
the /-mode perturbation field, which, in turn, is obtained 
from the decoupled field equations by standard 

numerical techniques. [The construction of F^^^ involves 
summation over the azimuthal number m, and — in the 
gravitational SF case — also over the ten tensorial har- 
monics for each multipole number I.] The direct piece 
F^" (and, hence, its I mode F^") is analyzed by local 
methods (see below). The expression (|l|) for the SF is 
then brought to the form [p|jl0| 



hm Fif ^Ao^L-B^-Ca/L 

x^zo 



-Do., (2) 



where hereafter the limit a; —> 2:0 is taken along the (in- 
going or outgoing) radial direction, L = l + 1/2, and A^, 
Ba, Con Etnd Da are certain "regularization parameters" 
(RP), derived from the local structure of F^^'- Roughly 
speaking, AaL + Ba + Ca/ L reflects the asymptotic form 
of F^Y [sind also F^"'') at large I, while Da is a residual 
quantity that arises in the summation over / (see [pUlO| 
for a more precise definition of the RP). Note that in the 
framework of the mode-sum method the task of calculat- 
ing the SF is divided into two separate parts: (i) Calcu- 
lating the quantities F^j^"^ — e.g., by numerically integrat- 
ing the decomposed field equations; and (ii) analytically 
calculating the four RP, by locally analyzing the direct 
part ijl^. It is the second part, the analytic derivation of 
the RP, that will primarily concern us in this Letter. 

In the original method by BO, the analysis of F^Y '^^^ 
carried out by applying a systematic perturbation expan- 
sion of the (Z-mode) Green's equation. The derivation of 
the full set of four RP required the calculation of the 
Green function up to second order in this expansion. So 
far, the RP where calculated using this original method 
in a few special cases: the scalar SF for a particle in radial 
or circular orbits and the gravitational SF for radial 
motion jl^,^ (all in Schwarzschild spacetime). How- 
ever, the required second-order perturbation analysis of 
the Green function turned out to be rather cumbersome, 
making it difficult to extend the calculation to more gen- 
eral cases. 

The goal of this Letter is to report on a new method 
for analytically deriving the four RP, which proves to be 
significantly simpler than the original BO's method. This 
new mode-sum method is based on an analytic expression 
constructed by Mino, Nakano and Sasaki [MNS] [^jl6| 
for the direct force F^" . (MNS developed this expression 
as part of a more comprehensive study aimed to provide a 
systematic description of the local structure of the metric 
perturbation.) Once this expression is available, essen- 
tially all we need is to decompose it into Z-modes, and 



then take the limit x ^ zq. The large-Z asymptotic be- 
havior of this expansion turns out to be a power series 
in 1/L, starting at L^, and its three leading-order coeffi- 
cients yield the three RP Aa , Ba , Ca according to 



lim F^Y = AaL + Ba + Ca/L + 0{L-^). 

x^zo 

The last parameter Da is then obtained by 
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Decomposing the right-hand side of Eq. (|^) and then sub- 
stituting Eq. (|^), one recovers Eq. (||). (Note that, since 
the tail part is regular, F^"" and F^"' admit the same 
large-/ asymptotic behavior, with the same coefHcients.) 

Applying this new method, we have been able to de- 
rive the values of all RP for an arbitrary geodesic orbit 
in Schwarzschild spacetime, for both the scalar and grav- 
itational self-force, as we now describe. Throughout this 
paper we use geometrized units G = c = 1 and metric 
signature — h-l— 1-, and denote the standard Schwarzschild 
coordinates by t, r, 9, ip. 

We begin by considering a simple scalar-field toy 
model. Once the scalar case is well understood, it be- 
comes fairly straightforward to obtain the RP in the grav- 
itational case as well. We shall consider a particle car- 
rying a scalar charge q, which (in the lack of SF) moves 
along an equatorial geodesic z^(t) of the Schwarzschild 
background (r being the proper time) with specific en- 
ergy and angular momentum £ and £, respectively. The 
scalar field <I> satisfies the field equation <I>.j^ = — 47rp, 
where p = qj 6'^{x — z){—g)~^/'^dT is the scalar charge 
density {g being the metric determinant). In this model, 
the force exerted on the scalar charge by an (external, 
regular) scalar field ^{x) is simply given by Fa = q^,a- 

Let e{x) denote the spatial distance from x^ to the 
worldline z^(t) (i.e., the length of the short geodesic con- 
necting to the worldline and normal to it). Using the 
Hadamard expansion of the field equation, MNS have 
been able to show that the "direct" scalar force takes the 
form 111 



F^ix) = q^ 



dir 



$^-(a;) ^ g//e. 



(5) 



where $'^"'(0::) is the so called "direct scalar field", and 
f{x) is some regular scalar function satisfying / = 1 + 
0{5x^), where 5x^ = x^ — Zq (the explicit form of / will 
not be needed here). For the analysis below it is useful 
to define S = e^. The direct force then takes the form 



F^'^\X) = q' [s-'^'la - fS-^/^S.a/i) 



(6) 



Being a regular function oi x, S may be expanded 
around zq as S* = 5*0 + 0{5x'^), where 5*0 = (.g^^ + 
u^u,^)Sx^Sx'^ , and u'^ = dz^ /dr. (In the expression for 
•^o, and Ua are to be evaluated at zq.) 
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In the next step we need to decompose the compo- 
nents F^" into spherical harmonics and then derive the 
four RP, as prescribed above. This goal was undertaken 
independently (and in two slightly different ways) by our 
two groups, BO and MNS, as we now describe. 

In the analysis by BO Q, the expression S = So + 
0{Sx^) is first substituted in Eq. With some simple 
manipulations, this equation can be brought to the form 

F^''i^) = ^o'Pi'^ + ^'P^ + ^'p2^ + 0{Sx), (7) 







Be) 



where eo = •S'o^^i ^-nd Pa^' denote terms of homogeneous 
order n in Sx". Notice that the three terms on the right 
hand side are of order 6x~'^, Sx~^, and Sx'^, respectively. 
The 0{Sx) correction term will not concern us here, as 
it does not contribute to the direct force at zq. 

As explained above, the mode I of F^" (or F^"'') is ob- 
tained by summing the contributions of all possible values 
of m for this specific /. This sum is invariant under rota- 
tion of the angular variables 6, ip. Using this invariance, 
BO chose a new, rotated, set of angular variables 9' , ip' 
(with respect to which the spherical harmonics are to be 
defined), such that the evaluation point zq is located at 
the pole, 9' = 0. Since all m 7^ spherical harmonics 
vanish at 9' = 0, with these rotated coordinates one only 
needs to consider the m — modes. This greatly simpli- 
fies the analysis. The presence of a coordinate singular- 
ity (the trivial 9' = singularity) at the evaluation point 
somewhat complicates the situation. To overcome this 
difficulty, BO introduced two regular "Cartesian-like" co- 
ordinates X, y. With this choice of coordinates {t, r, x, y), 
the /-decomposition of Eq. (0) becomes especially sim- 
ple. After taking the limit x — > zq, one finds that (i) 

f; (ii) the 



the term Pa^ yields zero contribution to F^ 



"^0 

term e^^P!^^ yields a constant (i.e., /-independent) con- 
tribution, which we denote 6q,; and (iii) the contribution 

of the term Pa is precisely proportional to L — we 
shall denote it OaL. From Eqs. (||) and (^) it now fol- 
lows that Aa — Oa, Ba = ba, and Ca — Da = 0. The 
detailed calculation of the coefRcients and ba is given 
in Ref. [Q. The final expressions obtained by BO for the 
RP (expressed in the original 9, (p coordinates, in which 
the motion is equatorial) are given by A'^g — Bg"^ — , 
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< = 0, 



Bl' 



(f_ {P - 2£^)k{w) + (r^ + £'^)E{w) 



g2 Sr^Kjw) - 2E{w)] 



g2 r\k{w) - E{w)\ 



^ r TT{C/r)V^n 



(8a) 
(8b) 
(8c) 
(8d) 



(with 'sc' signifying the scalar force RP), where K(w) 
and E{w) are the complete elliptic integrals of the first 
and second kinds, respectively, w = C? ji^C? -t- r^), / = 
(I-2M/r), y = l + £Vr2, andf = drjdT. The "±" sign 
in Aa refers to whether the limit x ^ zq \& taken along 
the ingoing or outgoing radial direction. The values (|^) 
agree with those obtained in j|] in the special cases L — ^ 
or f = 0. 

In an independent analysis, MNS directly used the 
"standard" 0, p coordinates (in which the motion is equa- 
torial) for decomposing the direct field. In this setup the 
TO 7^ modes contribute as well. MNS derived an an- 
alytic expression for the contribution of each /, m mode 
of the direct force, expanded in powers of Mjr. Then, 
MNS have shown the convergence of this expansion, and 
were able to explicitly sum it up (and sum over to), after 
which they recovered all RP values (||). 

We now turn to discuss the gravitational SF. Let 
K-^^ = ha0 denote the trace-reversed metric perturba- 
tion induced by a particle of mass /i (namely, haf3 = 
hai3 — \gap9^'^h^i,^ where hap is the metric perturbation 
itself and gap is the background metric). Just as in the 
scalar case, this "full" perturbation is the sum of the "di- 



rect" part and the "tail" part, /i^'^' 



h% + h'^l The 



gravitational SF is then obtained by applying a certain 
differential operator to h}ap'- 



TpscM T-itail — 7 /375 7 tail 



(9) 



evaluated at a; — > zq, where 



-u°'uf\'<u^ /2, (10) 



all quantities evaluated at zq. Let us also define 

p(dir,full) _ ,, i, /375r(dir,full) 



(11) 



after which Eq. ^ is recovered. 

The direct trace-reversed metric perturbation was 
found by MST . It can be expressed as 



h'p"{x) ^ A^.f{x,z)e ^ii/jUy, 



(12) 



where / is a regular function satisfying / = 1 -I- 0{Sx^) 
and iLa is the four- velocity parallelly-propagated from the 
particle's worldline to x along the short normal geodesic. 

As it stands, the tensor k is only defined on the par- 
ticle's worldline, because it involves u". In Eq. (|^) this 
on-worldline definition suffices, because h*^^^^ is regular 
at the particle's location However, /i^"'' and h'^" di- 
verge like 1/e (and their gradients like e~^). Therefore, 
when using Eq. ( pi] ) [e.g., when substituting it in Eq. 
(^] we must prescribe an extension of k off the world- 
line. Recall, however, that no ambiguity is caused by 
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this non-uniqueness of k: One just needs to use the same 
extension for both the direct and full forces. Below we 
consider several extensions: The one obtained by substi- 
tuting u" — > u" is denoted fc; and the extension based on 
"fixed (contravariant) components" (in the Schwarzschild 
coordinates) is denoted k. A third extension, mentioned 
below, is denoted k. Correspondingly, we shall denote 
the quantities associated with the above three extensions 
(e.g. the RP) by a hat, tilde, or an over-bar, respectively. 

The gravitational RP are obtained by decomposing the 
components F^" into scalar spherical harmonics, just as 
in the scalar case. In the fc-extension, MNS obtained 

= R^^ (with q ^ fj,), where Ra stands for all four RP, 
and the label "gr" signifies the gravitational case values. 
In the fc-extension, BO found @ Rg' = K^Rf, where 

For implementing the mode-sum method, Eq. (|^), we 
must have at hand i^^"", the spherical harmonic I mode 
of the (a component of the) full force. In the scalar case, 
this quantity is obtained by numerically computing the 
full-field mode <&; and differentiating it with respect to 
x" . In the gravitational case the procedure is a bit more 
involved. The linearized Einstein equation are separated 
by the decomposition 



(16) 
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hn i—1 



ir,t)Y^^'"\e,^), 



(13) 



where Y^^''™' are the ten (Zerilli-type) tensorial harmon- 



ics ||10[] . The quantities /i^')'™ can be computed by numer- 
ically integrating the decoupled field equations. Then, 



These values coincide, for C — 0, with the ones obtained 
]l2t using local analysis of the Green function. 

Let us now summarize the prescription for constructing 
the gravitational SF: (i) Numerically compute the radial 
functions /i*^*^' "^{r,t) (e.g., in the harmonic gauge Ipof); 
(ii) use Eq. (|l^) to construct the full modes F^^^; (iii) 
use Eq. ( p^ ) [along with Eq. (||), with q /i] to obtain 
the RP; and (iv) apply Eq. (|^). This prescription is now 
being implemented by Barack and Lousto . 

Finally, it should be reminded that the gravitational 
SF is a gauge-dependent notion, as discussed in Ref . . 
Nevertheless, the RP are gauge-independent The 
above prescription applies to the SF associated with the 
harmonic gauge, or any other gauge related to it by a 
regular gauge transformation. 
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To obtain Flf\ we need to decompose this quantity in 
ordinary (i.e., scalar) spherical harmonics Y^™{9, ip), and 
to sum over i and m. The outcome of this decomposition 
will depend on the extension of k. This decomposition 
appears to be difficult to implement with the fc-extension. 
For most extensions (including k) , infinite number of ten- 
sorial harmonics /i^*-'' will contribute to a single scalar- 
harmonic term F^^^ (due to the nontrivial dependence 
of k on 9). However, BO designed a third extension 
k = k + 5k, where 5k is a certain O(fe^) correction term 
(see ||l4j), in which only finite number {I — "i < V < I + 

of tensorial harmonics /i " contribute to a single term 
full j^^fuii givej^ by 



Ki'^n- Then 



pfull 
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i—1 V—l — 3ra— — l' 



where X*^*;^' ™ is a certain Ist-order differential operator 
(independent of 9, ip) given explicitly in Q. It can be 
shown ijl^ that the extension difference 5k, being 0{5x^), 
does not modify the RP values; hence. 
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